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Abstract. - We use computer simulations to study the microscopic dynamics of an athermal 
assembly of soft particles near the fluid-to-solid, jamming transition. Borrowing tools developed 
to study dynamic heterogeneity near glass transitions, we discover a number of original signatures 
of the jamming transition at the particle scale. We observe superdiffusive, spatially heterogeneous, 
and collective particle motion over a characteristic scale which displays a surprising non-monotonic 
behavior across the transition. In the solid phase, the dynamics is an intermittent succession of 
elastic deformations and plastic relaxations, which are both characterized by scale-free spatial 
correlations and system size dependent dynamic susceptibilities. Our results show that dynamic 
heterogeneities in dense athermal systems and glass-formers are very different, and shed light on 
recent experimental reports of 'anomalous' dynamical behavior near the jamming transition of 
granular and colloidal assemblies. 



Many materials, from emulsions and suspensions to 
foams and granular materials arc dense assemblies of non- 
Brownian particles pQ. Since thermal energy is irrelevant, 
the dynamics of these systems must be studied in driven 
non-equilibrium conditions. Together with the driving 
mechanism, a second important control parameter is the 
volume fraction of the particles, <p. These materials un- 
dergo a fluid-to-solid 'jamming' transition as <fi increases 
beyond some critical density (f) c [2]. 

The properties of this jamming transition have received 
considerable interest in recent years, and much progress 
was made through the analysis of idealized theoretical 
models, such as soft, frictionless, repulsive particles as 
studied below 3 . Structural and mechanical properties 
of systems on both sides of the jamming transition have 
been analyzed [4j[5] , and a number of remarkable features 
emerged, such as algebraic scaling of linear mechanical 
response [6], or the development of nontrivial timescales 
or lengthscales characterizing the macroscopic behavior of 
the system [7HS]- 

With the developments of new experimental techniques 
(confocal microscopy, original light scattering techniques, 



etc.) , the dynamics of systems near the jamming transition 
can now be resolved at the particle scale. Very recently, 
a number of 'anomalous' or 'unexpected' dynamic behav- 
iors were reported in granular and colloidal assemblies: 
superdiffusive particle motion [Tu]]T!l , nonmonotonic vari- 
ations of characteristic scales across the transition [TTHI4] , 
anomalous drop of dynamic correlations upon compres- 
sion [I3J[T5]- Thus, it appears timely to investigate, at 
the fundamental level of single particle trajectories, the 
signatures of the jamming transition in idealized theoret- 
ical models. Our work reveals original signatures of the 
jamming transition in both fluid and solid phases which, 
we believe, shed light on recent experimental findings. We 
also establish a connection between single particle dynam- 
ics and collective particle motion, which allows us to de- 
velop an intuitive and appealing picture of jamming as the 
consequence of the diverging size of rigid particle clusters. 

We study a bidimensional 50:50 binary mixture of 
iV particles of diameter ratio 1.4 with harmonic repul- 
sion [T2], Vljij) — k(rij — CFij) 2 , where is the distance 
between the centers of particle i and j, Cij = (oi + aj)/2, 
and Ci is the diameter of particle i. Particles only inter- 
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act when they overlap, V(rij > <7ij) = 0. We work in 
the zero temperature limit, so that all the dynamical pro- 
cesses in the system are induced by an external driving. 
We use a simple shear flow, and perform quasi-static shear 
simulations, as described before [T7]- Very small shear 
strains are applied, followed by a minimization of the po- 
tential energy. At each step, particles are first affinely 
displaced along the y-axis by an amount Sy = jqx with 
70 = 5 x 10 -5 . Appropriate Lees-Edwards periodic bound- 
ary conditions are used. During the subsequent energy 
minimization, additional nonaffine displacements occur in 
both directions. In the dilute limit, these nonaffine dis- 
placements are absent, and their presence directly reveals 
the influence of interactions and interparticle correlations. 
In the following, we focus on the purely nonaffine displace- 
ments occurring along the x-direction, i.e. transverse to 
the flow. We use system sizes N = 900, 1600 and 2500 to 
detect finite size effects. The unique control parameter is 
the volume fraction, <j>. Below a critical value, (j> c ~ 0.843, 
which is well defined and sharp in the thermodynamic 
limit, the system flows with no resistance, while a yield 
stress grows continuously from 0, when <j) increases above 
4>c |18j . This corresponds to the jamming transition for 
the present system and driving conditions. 

The simplest quantity which is measured from parti- 
cle displacements is the root mean-squared displacement, 

1 /2 

A(7) = (w E,=iM7) - ^(0)] 2 ) , where a* (7) - a* (0) 
is the transverse displacement of particle i after a total 
strain 7 is applied. In this expression, the brackets in- 
dicate averages taken over independent initial conditions 
all chosen in the driven stationary state at a particular 
density. We present the evolution of A (7) at various vol- 
ume fractions across <j) c in Fig. [T] The displacements in 
the yield-stress flow regime above </> c are diffusive at all 
strains, A ~ ^7. This is a robust signal of the quasi- 
static dynamics in yield-stress materials. It can be traced 
back to the accumulation of plastic rearrangements span- 
ning the entire system [191120) . The important observation 
here is that the diffusion constant depends only weakly on 
density and distance to <f> c . 

The regime below 4> c is more interesting. At any given 
volume fraction, particles move in a superdiffusive, ballis- 
tic manner at short strain, A ~ 7, before crossing over 
to diffusive motion at large strain. Since we work at zero 
temperature, this ballistic regime is completely unrelated 
to the trivial short-time ballistic displacements in particle 
systems with classical Newtonian dynamics. 

Superdiffusive behavior can be interpreted from the ob- 
servation that below the jamming transition, the network 
of particle contacts is not sufficient to insure mechan- 
ical stability, and there exists a large number of zero- 
frequency modes allowing particle displacements at no en- 
ergy cost [21] . The actual displacement field is a superpo- 
sition of these modes, which can thus persist until the par- 
ticle configuration has evolved sufficiently to decorrelate 
the mode spectrum. Persistence of the modes, followed 
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Fig. 1: Nonaffine displacement A (7) for various volume frac- 
tions (f>. For <j> > 4> c , diffusive motion is observed at all 7. In 
contrast, a superdiffusive regime develops for <f) < <j> c , with a 
particle 'velocity' £a(4>) = A(t)/7 which increases when ap- 
proaching <f> c . 

by their decorrelation explain superdiffusive and diffusive 
motion, respectively. 

Remarkably, Fig. []] shows that particles move faster 
when is increased towards <f> c in both ballistic and diffu- 
sive regimes. In particular, the short-time particle 'veloc- 
ities' £a(<P) = hm 7 _j.o A(7)/7 increase with 4>. This is a 
surprising finding because compression towards jamming 
usually yields slower dynamics [TTJ[TS]. Note that I a is a 
displacement per unit of strain and thus has the dimen- 
sion of length. Thus an increasing particle velocity in fact 
implies an increasing length scale. Indeed, we will argue 
below that superdiffusion reflects the concerted motion of 
solid-like clusters of particles correlated over a length scale 
comparable to ^a(<^) ; and which grows towards <p c . 

The displacements A and the associated length scale 
£a{4>) measure how much particles have to move in addi- 
tion to the affine displacements, in order to accomodate 
the imposed shear flow. At cf> = 0.840, for example, this 
additional displacement is about ten times larger than the 
affine contribution itself. This means, that close to <fi c 
the system is in a highly 'fragile' state and small changes 
in the boundary or loading conditions lead to large-scale 
motions [22] . 

Interestingly, superdiffusion was recently used to iden- 
tify the location of the jamming transition in a bidimen- 
sional granular system [11]. Due to the 'random agitation' 
driving mechanism in that experiment, superdiffusion is 
observed on a modest time window, and only very close to 
4> c . This is consistent with our finding that superdiffusion 
is indeed most pronounced near the transition. 

A second, often studied, correlation function to quantify 
single particle dynamics is the 'overlap' function [TT] 

«^ = (lf> P (-M2>_^. (1) 
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Fig. 2: The average overlap Q(a,y) for a — O.OOf and different 
volume fractions across the jamming transition. The relax- 
ation becomes faster when <j> increases towards <f> c , but slows 
down when <j> increases above (j> c . The fluid and solid phases 
are thus readily distinguished by their distinct behavior upon 
compression. 

The overlap Q(a, 7) goes from 1 to as typical particle dis- 
placements get larger than the probing length scale a. It is 
thus very similar to a self-intermediate scattering function 
in liquids, with a playing the role of an inverse wavevector, 
and 7 the role of time. 

The behavior of Q(a, 7) is presented in Fig.[2]for a fixed 
a = 0.001 (qualitatively similar results are obtained for 
different choices of a) and various volume fractions across 
4> c . For (f> < 4> Cl we find that the overlap function de- 
cays faster when <p increases. This is consistent with the 
above observation that displacements get larger closer to 
4> c - A qualitatively different behaviour is found above the 
transition, where the overlap decays more slowly when <p 
increases. This is in striking contrast with the behavior of 
the mean squared displacement in Fig. [T] which showed no 
such variations with <j). 

The data in Fig. [2J thus imply the existence of an un- 
expected nonmonotonic variation of the dynamical be- 
haviour, which is absent in the root mean-squared dis- 
placement, Fig. [TJ To quantify these differences we define 
the analog £q(4>) of the displacement scale, £&{(j)), dis- 
cussed above. For a given strain 7 we measure the value 
of a such that 

Q (a = 7^(0), 7) =0.5. (2) 

The value 0.5 is arbitrary and we find similar results with 
0.3 and 0.7. For small strains within the ballistic regime 
this definition insures that £q is independent of strain 7 
and reflects a short-time velocity, as does its analog £a- 
For larger strains a crossover to diffusive behavior sets 
in similar as for the root mean-square displacement A 
(Fig- HI- In the following, we want to avoid this crossover 
and thus choose 7 as small as possible. 

In Fig. [3] we report the behaviour of £q{4>) as denned 



Fig. 3: Relaxation length scale £q(4>) quantifying the decay 
of the overlap Q(a,'y) for a single strain-step 7 = 70, from 
Eq. ([2]) (£q is defined from purely elastic deformations above 
jamming). It is a nonmonotonic function, with a maximum at 
the jamming transition, mirroring the behavior of the overlap 
in Fig. [2] It obeys power law behavior with different exponents 
on both sides of the transition. In the fluid phase, £q behaves 
similarly to the 'velocity' £a defined from superdiffusive be- 
haviour. 

from Eq. ([2]), which clearly reflects the nonmonotonic be- 
haviour of the overlap observed in Fig. [2] This result shows 
that a remarkably simple statistical analysis of single par- 
ticle displacements in athermal systems very easily reveals 
the existence of, and quite accurately locates, the fluid-to- 
solid jamming transition. This measurement requires no 
finite size scaling, or other involved or indirect statistical 
analysis [HMOS]. 

We have discussed two signatures of the jamming tran- 
sition from the behavior of A (superdiffusion) and Q(a, 7) 
(nonmonotonic relaxation as a function of density). It 
should come as a third surprise that both dynamic quan- 
tities apparently provide distinct informations, in partic- 
ular for <p > <j) c , although they both aim at quantifying 
dynamics at the particle scale. In fact, such discrepancies 
are well-known and well-studied in the field of the glass 
transition where several similar 'decoupling' phenomena 
have been studied [23] • Indeed, if the particle dynamics 
were a Gaussian process, the information content of A 
and Q would be mathematically equivalent. Their differ- 
ent behavior thus suggests that the distribution of single 
particle displacements is strongly non-Gaussian [25] . 

To substantiate this claim we report the evo- 
lution of the van Hove distribution, P(x,^/) = 
{jt $ (x — [xi (7) — Xj(0)])), as a function of volume 
fraction below <fi c (Fig. [4j . Above <p c the shape of the van 
Hove function has a form well-known from recent work 
on sheared glasses [H2[20]; therefore, we do not show our 
data here. Briefly, the distribution is 'bimodal', with a 
first component at small x resulting from displacements 
during reversible elastic branches, and a second compo- 
nent at larger x due to irreversible plastic events. This 
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Fig. 4: Displacement distribution after a single strain step, 7 = 
70, for different volume fractions </> < <j>c, rescaled such that a 
Gaussian process (dashed line) has unit variance. Increasingly 
heterogeneous distributions of particle displacements are found 
when approaching <f) c , with the development of algebraic tails 
yielding a diverging kurtosis at <j> c . 

latter contribution completely dominates the second mo- 
ment of the distribution, and thus controls the diffusive 
behavior of A (7) in Fig. [1] 

The decomposition between elastic and plastic branches 
above <j> c suggests a separate analysis of these two types 
of dynamics. We have computed Q(a, r y) along elastic 
branches only, and defined £q (</>), the 'elastic' analog of 
the displacement scale £q(4>) for the full overlap; it is 
shown in Fig. [3] Both £q and £q are very close, con- 
firming that Q(a, 7) is not very sensitive to plastic rear- 
rangements, in contrast to A(7). The scaling of the typi- 
cal elastic response above 4> c has recently been discussed 
in great detail [24l [26l [27] , and we simply quote the re- 
sult: £q ~ ((j) c — (j))^ 1 ^ ~ £q, in good agreement with our 
numerical results. 

Below (p c the distributions remain highly non-Gaussian 
(Fig.|4|), even though plastic events are absent. Deviations 
from a Gaussian shape become stronger as <j) increases to- 
wards 4> c . The distribution develops polynomial tails, and 
it is thus not surprising that different averages taken over 
such distributions provide quantitatively distinct results. 
In Fig. El we show that t Q ~ {<j> c - 0)~ a9 , for 4> < <j> c . 
In contrast, we find that the scale £a, defined in Fig. [IJ 
scales as £a ~ (</> c — 4>)~ XA with an exponent slightly dis- 
tinct from the one of £q. Hence, despite the appeal of the 
algebraic 'scaling' relations reported in Fig. [3l our parti- 
cle distributions P(x, 7) do not follow simple scaling forms, 
even along elastic branches, in apparent contradiction with 
recent claims |27j . Single particle dynamics is not 'univer- 
sal' near the jamming transition and it is somewhat am- 
biguous to define 'typical' displacements from such broad 
distributions. 

It is tempting to speculate that the polynomial tails 
which develop near but below <p c in Fig. [4] are in fact the 
'precursors' of plastic events taking place above <j) c . The 




Fig. 5: Grey-scale plot of the spatial fluctuations of the overlap 
in Eq. J]} (dark='immobile', light='mobile') for (f> = 0.825 < 
<j) c (left), (j} c = 0.840 (middle) and <f> = 0.85 > <j> c (right). Strain 
7 = 70 and probing length scale a chosen such that Q « 0.5. 
Spatial correlations are increasing with cj> towards <f> c , but stay 
comparable to system size above <f> c when the system is in the 
solid phase. 

tails in fact become so broad that the kurtosis of the distri- 
bution (also called 'non-Gaussian parameter') actually di- 
verges at small 7 when </> c is approached from below (data 
not shown). This implies that the kurtosis could be used 
as a second, simple statistical indicator of the underlying 
jamming transition. In glass-formers, the non-Gaussian 
parameter increases as the glass transition is approached, 
but much more modestly |28j , while van hove distributions 
develop exponential rather than algebraic tails |25j . 

Let us now turn to the discussion of particle correla- 
tions. This will allow us to make concrete the notion of 
'rigid' clusters that we have alluded to above. To this end, 
we first discuss the images shown in Fig. [SJ where the spa- 
tial fluctuations of the overlap are shown. From this figure, 
we can readily identify dynamically correlated clusters of 
particles, that clearly grow in size upon approaching </> c 
from below. These images are direct evidence that dynam- 
ics becomes more collective as <j) c is approached. Above </> c 
the system is now solid, and responds as an elastic body. 
Correspondingly, mobility fluctuations are correlated on a 
length scale comparable to the system size, the only visi- 
ble effect of the density being that the snapshot looks less 
'disordered' at the largest <f>. Very similar observations 
were recently made experimentally in a soft granular ma- 
terial [T4] . 

To quantify these qualitative observations, we study the 
dynamical susceptibility Xii defined as the variance of 
statistical fluctuations of the overlap Q X4( a i7) = 
N ^Q(a, 7)^ — Q(a, 7)^ , where Q represents the in- 
stantaneous contribution to the average Q. Our results 
are summarized in Fig. [6j In all cases, Xi( a >j) displays 
a well-defined maximum when shown as a function of 7, 
which shifts towards lower strain when <j) increases. This 
maximum identifies a strain for which fluctuations of the 
overlap are maximal, and is slaved to the typical strain 
over which Q decays in Fig. [21 as commonly found in stud- 
ies of dynamic heterogeneity [30|l3Tj . As is well-known, the 
height of this maximum is a good estimate of the number 
of particles relaxing in a correlated manner |31j . which jus- 
tifies the fundamental importance of xt- The increase of 
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Fig. 6: Dynamical susceptibility x& as function of strain 7 for: 
(Top) different volume-fractions cj> < <j> c and probing length 
scale a. (Bottom) a single <j> — 0.85 > <f> c and different system 
sizes. The dynamic susceptibility increases when <j> increases 
towards <f> c reflecting increasingly collective dynamics. Above 
cf> c dynamic correlations are system size dependent, even when 
only elastic branches are considered. 



the height of the peak of the dynamic susceptibility when 
4> increases towards 4> c in Fig. [5] thus quantitatively con- 
firms the visual impression given by the snapshots: the 
dynamics is increasingly collective when increases to- 
wards <j) c . As shown in Fig. [SJ this conclusion does not 
depend on our choice of a probing length scale a. 

To quantify the possible divergence of the spatial corre- 
lation of the dynamics, we analyze the data as above for 
the overlap. For a given strain 7, we look for a such that 
X4 is maximum (it corresponds roughly to £q( , j) defined 
in Eq. ©) and measure its height. We report the evo- 
lution of this maximum as a function of volume fraction 
in Fig. [71 The data for three different system sizes super- 
impose below (f> c , and we can satisfactorily describe the 
^-dependence as: xa ~ {4>c — </>)~ 1,8 - Comparing this ap- 
parent divergence with the much more modest increase of 
dynamic correlations in glass-forming liquids close to the 
glass transition [32] , emphasizes once more the qualitative 
differences between both types of transitions. 
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Fig. 7: Evolution of the maximum of the dynamic susceptibility 
X4 across the jamming transition (7 = 70; symbols as in Fig. [6j . 
It diverges algebraically as <j> increases towards <j>c and is system 
size independent. Above cj> c , the maximum is proportional to 
the system size, with an amplitude which decays slowly upon 
compression, both for plastic and elastic dynamics. 



Assuming that correlated domains have a compact ge- 
ometry, this finding immediately translates into a gen- 
uine dynamic correlation length scale, £4, which grows as 
£4 = \/X4 ~ (<Pc ~ (f>)~°' 9 - Strikingly, we found a similar 
power law divergence for both £q and £\ in Figs. [T] and [H 
The exponents are furthermore close to the one obtained 
in previous work [17] using the spatial dependence of the 
displacement correlation function. 

These findings call for a simple relation between single- 
particle displacements, as characterized by £q or £/\, and 
particle correlations, as given by £4 or by displacement cor- 
relations. The picture that we propose assumes particles 
to form 'temporarily rigid' clusters of the size of the corre- 
lation length £4. Driven by the shear-strain 7 these clus- 
ters move or rotate in a solid-like manner over distances 
A c i us tcr = £4 - 7- This motion shows up as pronounced bal- 
listic regime in the mean-squared displacement or in the 
decay of the overlap function Q and allows to make the 
identification, £q ~ £4. The cross-over to particle diffusion 
then corresponds to a typical cluster lifetime. On longer 
strains, clusters break up and lose their identity. This 
process is evident in the reduction of the dynamical cor- 
relation length as measured by the decreasing amplitude 
of Xi with 7 in Fig. [6] 

Moving to <f) > 4> c in Fig. [6l we observe that \4 still ex- 
hibits a maximum at a given strain 7, which again tracks 
the relaxation of Q. Remarkably, we now find that the 
height of this peak strongly depends on the system size, 
and increases roughly linearly with N, see Fig. [6] This 
suggests that spatial correlations saturate to the system 
size at 4> c and remain scale-free above <p c in the entire 
solid phase [T5]. This is in agreement with studies of the 
mechanical response of disordered solids, which showed ev- 
idence of plastic rearrangements spanning the entire sys- 
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tern [2U][55]- Additionally, we show in Figs. [5] and [7] that 
the dynamic susceptibility measured along purely elastic 
branches (discarding plastic events), is linearly growing 
with N as well. Indeed, the elastic response of a disor- 
dered solid is also scale- free, as shown theoretically [54 , 35 
and numerically [551157] . 

For a given system size, %4 grows with <j> below <p c , and 
becomes proportional to N above <j> c . However, the am- 
plitude of Xi decreases slowly with <f) above </> c . The net 
result is that %4, for a fixed system size N, has a striking 
nonmonotonic behaviour with <fi, and presents an abso- 
lute maximum at 4>c > which adds to our list of microscopic 
signatures of the jamming transition. 

A similar maximum of Xi nas recently been reported 
both for colloidal [T51IT3] and granular [TT] assemblies, and 
given two distinct interpretations. It was associated to a 
nonmonotonic dynamic correlation length scale £4 diverg- 
ing at 4> c in Ref. [TT] . Alternatively, it was attributed to a 
nonmonotonic strength of essentially scale-free spatial cor- 
relations in Refs. [T3lfT5] . Our results clearly favor the sec- 
ond interpretation to explain the behavior above <p c . The 
dynamic susceptibility Xi above 4> c is system size depen- 
dent with scale- free spatial dynamic fluctuations (Fig. [TJ. 

This finding suggests that the length scale identified 
from the low-frequency part of the vibrational spectrum, 
and much studied recently [7]I27II37| does not directly influ- 
ence the dynamics at the particle scale beyond the linear 
regime where it belongs. We suggest, however, that it in- 
directly influences the behaviour of xa, and leads to the 
decaying strength of the correlations as observed in Fig. [7] 
above <p c . The evolution of the displacement length scale 
£q above <f> c in Fig. [3] suggests a simple explanation for 
this decay. The overlap function Q(a, 7) relax when typi- 
cal displacements are of the order of a. When <f> increases, 
typical displacements in response to strain increments be- 
come smaller, and more steps are needed to decorrelate the 
overlap. If these steps are not perfectly correlated, more 
steps directly imply less dynamic fluctuations and thus a 
reduced Xi [31] . This physical explanation was made more 
quantitative in the empirical model of Ref. [15] , which then 
accounts for a nonmonotonic evolution of Xi across <p c . 

To summarize, borrowing tools from studies of dynamic 
heterogeneity in glass-formers, we provided a detailed ac- 
count of the evolution of the microscopic dynamics of an 
idealized model system across the jamming transition. We 
found a number of original signatures, that are strikingly 
different from the behaviour of viscous liquids, but might 
have been observed in recent experimental reports. Given 
the (relative) ease with which particle trajectories can 
be monitored in suspensions, foams, or granular assem- 
blies, we hope our study will motivate further experimen- 
tal studies along the lines of the present study. 
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